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We revisit the classical but as yet unresolved problem of predicting the breaking on¬ 
set of 2D and 3D irrotational gravity water waves. This study focuses on domains with 
flat bottom topography and conditions ranging from deep to intermediate depth (depth 
to wavelength ratio from 1 to 0.2). Our calculations based on a fully nonlinear bound¬ 
ary element model investigated geometric, kinematic and energetic differences between 
maximally recurrent and marginally breaking waves in focusing wave groups. Maximally 
steep non-breaking (maximally recurrent) waves are clearly separated from marginally 
breaking waves by their normalised energy fluxes localized near the crest region. On the 
surface, this reduces to the local ratio of the energy flux velocity (here the fluid velocity) 
to the crest point velocity for the tallest wave in the evolving group. This provides a 
robust threshold parameter for breaking onset for 2D and 3D wave packets propagating 
in uniform water depths from deep to intermediate. Warning of imminent breaking onset 
was found to be detected up to a fifth of a carrier wave period prior to a breaking event. 

Key words: Authors should not enter keywords on the manuscript, as these must be cho¬ 
sen by the author during the online submission process and will then be added during the 
typesetting process (see http://journals.Cambridge.org/data/relatedlink/jfm-keywords.pdf 
for the full list) 


1. Introduction 

Despite its long research history, the physics underpinning the breaking of water waves 
has remained incompletely understood, including prediction of its onset and strength. 
Yet this knowledge is of fundamental importance in quantifying atmosphere-ocean ex¬ 
changes, determining structural loadings on ships and platforms, and optimising opera¬ 
tional strategies for maritime enterprises. 

Many criteria for pre dicting b reakin g onset of water waves have been proposed since 
the pioneering study of StokesI ( 1847f) . These criteria arise from theoretical arguments 
based on idealized models, numerical simulations, laboratory experiments and field ob¬ 
servations. However, while adding many insights, these approaches have not yielded a 
robust breaking threshold for phase-resolved waves in the physical domain, reflecting the 


f Email address for correspondence: x.barthelemy@unsw.edu.au 










2 


X. Barthelemy, M.L. Banner, W.L. Peirson, F. Fedele, M. Allis, F. Dias 


complexity of the underlying dynamical processes. In fact, there is a glaring absence of 
a precise definition of breaking onset. 


Briefly, it has long been considered that breaking is a process with a threshold, with 
criteria for predicting breaking onset falling into three categories: geometric, kinematic 
and energetic. The majority of breaking criteria have been based on a geometric or kine¬ 
matic threshold, and mainly limited to plane (2D) waves. Geometric threshold variables 
have included wave steepness, wave asymmetry, maximum theoretical (global) steep¬ 
ness and the occurrence of a transient vertical segment on the forward face of the wave 
crest; kinematic threshold variables have included the Lagrangian crest acceleration and 
the ratio between c rest fluid speed and phase speed. The recent comprehensive review by 
Berlin et al. ( 201,11 1 provides an excellent overview of the collective observational and the¬ 


oretical effort and ou t comes based on kinematic/ g eomet r ic approaches. The rece n t con- 
tributions of Sherner (l2013h . Shemer fc Liberzon ( 2014 1. Kurnia fc van Groesen ( 2014 1 
and IShemer fc Ee ( 2015 1 add to this otherwise exhaustive coverage. Overall, current 
knowledge does not support a kinematic or geometric criterion that provides a generic 
threshold which differentiates breaking from recurrent behaviour for deep water waves. 
While the vertical tangent segment and kinematic criteria provide valid a posteriori con¬ 
ditions for breaking onset, they provide no dynamical insight or advance warning of 
imminent breaking. 


A third approach based on dynamical criteria has been explored to explain the onset 
of breaking. This concept is based on the evolution of the intra-group energy flux which 
causes the tallest crest of an unsteady wave group to break when a local stability threshold 
is exceeded. Monitoring the energy flux field in this highly nonline ar, unsteady flow 
enviro nment makes rigorous analysis difficult. The overview article bv lTulin fc Landrini 
( 2000ll highlights the very insightful inroads made by Tulin and his collaborators over 
the previous decade into unsteady nonlinear wave group evolution and breaking, based 
on intra-group energy flux theory, simulations, observations and analyses. One of the key 
results they proposed from their studies is that breaking onset is initiated within a wave 
group when a crest particle speed exceeds the linear group speed. Pending verification of 
its general validity, this criterion is able to signal breaking onset much earlier than the 
traditional kinematic criterion. 


Subse quently, BannerfcTian ( 1998h . ISong fc Banner! (2002) and the experimental 
study of Banner fc PeirsonI (j2007l l investigated a growth rate based on a parametric en¬ 
ergy convergence rate for 2D wave gro ups, using a frame of reference that tracks the wave 


group maximum. Berlin et al\ (l2013l) discuss the merits of this appr oach based on the 


furthe r study of Tian et al. (2008) for 2D wave breaking. Very recently, Derakhti fc Kirb^ 
( 201(ih reported very encouraging support for this approach in their numerical study of 
unsteady 2D wave packets in a model framework that can accommodate sequential (mul¬ 
tiple) breaking events as the packet evolves. They confirm ed the presenc e of sy stematic 
crest/trough leaning motions, as investigated in detail in Banner et al. ( 20141) . In the 
present study, we revisit this local energy growth rate approach for both 2D and 3D 
breaking onset simulations, for which our findings are given below in subsection 14.31 


Finally, significant additional challenges arise in representing wave breaking in broad- 
banded directional sea states in the spectral (wavenumber- frequency) context. This is t he 
domain used for computing ocean wave forecasts (e.g. see iGhalikov BabaninI ( 20121) ). 
In that context, there is even less consensus on how to predict/identify breaking events. 
However, this is beyond the scope of the present paper, where we focus on wave breaking 
in the physical space-time domain. 
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2. Rationale underpinning our breaking onset threshold investigation 

Present understanding of the physics underpinning wave breaking onset in the physical 
domain is fragmentary, including a precise definition. This has precluded reliable predic¬ 
tion of wave breaking onset even in controlled laboratory and numerical wave basin 
conditions. Our investigation directly addresses this time-honoured knowledge gap. 

Based on energy flux considerations, we propose a new breaking onset threshold pa¬ 
rameter. The behaviour of this parameter in the wave crest region is found to be of central 
importance. A new definition for a generic breaking onset threshold emerges naturally. 
Through an ensemble of simulations of diverse nonlinear wave packets in both deep and 
intermediate depth water over a flat bottom, we establish the existence of a breaking 
onset threshold band for this parameter, determining its upper and lower bounds. Below 
the lower bound of the proposed breaking threshold band, steepening carrier waves evolve 
through the packet envelope maximum without the occurrence of a vertical tangent in 
the wave surface profile. All maximally steep non-breaking carrier waves exist below this 
lower bound, which we also refer to as the maximum recurrence threshold. 

Thereafter, the smallest increment in the crest wave energy density (e.g., as produced 
by increasing the wave paddle amplitude) causes the breaking threshold parameter at 
the tallest crest to increase. After it evolves further, it exceeds our proposed breaking 



actively breaking crest depends on the wave scale, the influence of surface tension and 
on the strength of the breaking, which reflects the energy convergence rate at the wave 
crest. The strength of the breaking event is an unknown function of the magnitude of the 
breaking parameter above the threshold, with the smallest exceedance margins associated 
with very weak (i.e. marginal) breaking. 

For long carrier wavelengths, the evolution to breaking leads to plunging jets of varying 
size relative to the wavelength: the lower the energy input rate to the wave, the smaller 
the plunging jet. For short wavelengths, the jet formation is modified by surface tension 
and has been in vestigated in detai l with boundary element calculations, theory, and 
experiments (e.g. Tulin fc Landrinil (200^). When the wavelength exceeds about 2 m, 
breaking starts with the formation of a small plunging jet, just as it does when the 
surface tension vanishes. As the wavelength is decreased, surface tension forces become 
relatively larger and the jet tip becomes rounded. For wavelengths less than about 0.5 m, 
the jet is replaced by a bulge, and capillary waves appear upstream of the leading edge 
(toe) of the bulge, as shown in Figure 1 of Duncan (2001). In the present investigation, 
surface tension is not included explicitly. However, we show below (in section 15.51) that 
its estimated effect on our proposed generic breaking onset threshold is negligible for 
wavelengths longer than 1 m. 

In this study, we used a fully nonlinear wave code capable of capturing the initial stage 
of crest overturning, including the critical visible signature of a transient vertical tangent 
on the forward face of the wave crest. The results from this code were used to determine 
the upper and lower bounds of our proposed breaking onset threshold band. Note that 
we do not directly solve the instability problem for breaking onset, nor does our model 
provide information beyond the initial stage of crest overturning for waves that exceed our 
breaking threshold. In common with many other studies, this choice of using strongly 
chirped wave packets which focus rapidly minimises the evolution time and fetch to 
breaking, or recurrence, and hence the CPU time. It also reduces the adverse Lagrangian 
drift implications for the computational grid (see end of sectionfor details). However, 
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the high chirp rate also restricted our attention to weak breaking cases, as further small 
increments in paddle amplitude only produced wave breaking at the paddle. 

Our breaking onset parameter is formulated as the local energy flux relative to the 
local energy density, normalised by the local crest speed, and is operative throughout the 
subsurface region including the wave surface. For the condition of zero surface pressure, 
its projection on to the wave surface reduces to a simple quasi-kinematic form involving 
the ratio of the surface fluid speed u to the crest point speed c. We note that the bound¬ 
ary geometry does not enter the breaking onset criterion explicitly, so it is potentially 
applicable to variable depth bottom topography scenarios. 

In the present paper, the behaviour of our breaking onset threshold parameter is de¬ 
termined from an ensemble of numerical simulations of 2D and 3D chirped focusing wave 
packets on deep and intermediate depth flat bottom topography, ranging down to one 
fifth of the dominant wavelength. For this class of wave packets, the results establish the 
existence of a generic narrow threshold band for breaking onset, for which u/c is found 
to be appreciably lower than the traditional kinematic breaking criterion of u/c > 1. 
As foreshadowed above, the resemblance of the surface projection of our breaking onset 
threshold to a kinematic criterion is fortuitous. It only takes this form for zero surface 
pressure forcing conditions. Its intrinsic dynamical nature is confirmed by two additional 
factors: the concomitant subsurface threshold does not reduce to a kinematic form and 
the breaking threshold u/c ratio is considerably below unity. 

The findings from our numerical investigation for 2D deep and intermediate depth 
water waves are f ound to agree closely wit h measurements from our companion observa¬ 
tional studies by Saket etall ( 2017ah and Saket et al. ( 2017&I1 . Their observations also 
validated our proposed breaking onset threshold for moderate wind forcing and also for 
modulationally-focusing bimodal wave packets. 


3. Methodology 

3.1. Wave generation 

Wave groups, either 2D or 3D and in deep or intermediate depth water over flat bottom 
topography, were generated using a bottom-hinged flap-type snake wavemaker at one end 
of the tank. The motion Xp{t, y) of the wavemaker flap at the lateral location y followed 
the Class 3 ’chirp packet’ motion from ISong fc Banner! ( 2002l l (this is implicit hereafter 
in the C3 designator in each documented run hie): 


Xp(t,y) = -0.25Apl 1 


tanh 
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where t is the time, Ap is the amplitude of the paddle motion, N is the number of waves 
in the temporal wave packet, Wp is the baseline driving frequency of the paddle, with 
corresponding linear wavenumber fcp, and Cch = 1.0112 x 10“^ specihes the chirp rate 
used in this study. The phase d> (y, Xrrm,,,Y r r,nv) s p ecihes the coordin ates of the point of 
linear convergence Isee iDalrvmple fc Kirbvl (llQSSh : IPalrymplj ( 1989lD : 
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where 9 is the focal angle at location y along the paddle. The down stream boundary 
oppos ite the wave paddle is a fully-absorbing boundary condition, as in Grilli fc Horrillol 
( 199711 . 


3.2. Numerical wave tank 

There has been growing interest in the development of three-dimensional models which in¬ 
herently incorporate nonlinearity and associated dispersion effects. The broad-bandedness 
in both fre quency and direction o f real sea states poses significant challenges in numerical 


m both tre quency and direction o t real sea states poses sigmhcant challenges m numerical 
simulation. iBateman et~^. I (l200lh demonstrated the importance of directionality and the 


consequent benefits of efficient wa ve modelling when compari ng numerical simulations 
with the laboratory observations of .Tohannessen fc Swan ( 200lh . High-order spectral ex¬ 
pansion approa ches using effici e nt FF T sol vers for appl i cation to 3D waves have been 


developed (e.g. iDucrozet et al\ (|2012l l and iFedele et al. 


solve the full Navier-Stokes equations (jPark et al\ (|2003M . but viscous flow solvers tend 


( 201(ih l. Another option is to 


to be too dissipative and computationally time-consuming. 

Numerical models of 3D potential flow wave propagation can b e divided into thre e 
main categories : (a) boundary element integral methods (BEM): e .g. Baker_et_a^ 1982), 
Balemnjj_et_a^ 20^ 1. Clamond fc Grue JboOlh. Grilli et aL (I200D . Fochesa to et al. ( 2007 1. 


Guvenne &: Grilli ( 200dl . Xue et al.l ~ 2001 1. Hou fc Zhand ( 200^ Fructus et ai ( 2005 1: 


(b) fin i te element method (FEM) e.g. M^^etgL ( 20011; (c) spectra l methods: e.g. Dommermuth fc Yud 
( 1987l l. West et al. ( 1987 1. Craig fc SulemI ( 1993 1. Nicholy ( 1998ll . Bateman et al. ( 2001 1. 


Spectral methods based on perturbation expansions are known to be very efficient. 
These methods reduce the water-wave problem from one posed inside the entire fluid 
domain to one posed on the boundary alone, thus reducing the dimension of the formula¬ 
tion. This reduction can be accomplished by using integral equations over the boundary of 
the domain (so-called boundary integral methods) or by introducing bou ndary quantities 
which can be expanded as Taylor series for reference domai n geometrie s ( Xu fc Giivennel 


( 2 OO 9 I II. Both approaches have been summarised recently bv iMal ( 201(tl . BEM techniques 


are efficie nt for representing wave prop agation and overturning until the wave surface re¬ 
connects ( Grilli fc Subramanval ( 1996 


The present study used a boundary element numerical wave tank (he reafter NWT) code 
called WSIM, which is a 3D extension of the 2D code developed bv iGrilli et all ( 1989[l 


to solve the single-phase wave motions of a perfect fluid. It has been applied extensively 
to the solut ion of finit e amplitude wave propagation and wave breaking problems (see 


chapter 3 of iMal ( 2010ll l. 


The perfect fluid assumption makes WSIM unable to simulate breaking impact sub¬ 
sequent to surface reconnection. However, its potential theory formulation enables it to 
simulate wave propagation in a CPU-efficient way, without the diffusion issues of viscous 
numerical codes. The simulation of wave generation and developme nt of the onset of 
break in g events can be carried out with great precision, as shown by Eochesato fc Diai 
( 2006ll . Eochesato et al. ( 2007h . 


depth water and shows excellent energy conservation (Grilli et al. 

(1989ll. Grilli & Svendsen 

(Il99f)ll. Grilli & Subramanva |l994ll. Grilli & Subramanval (1 19961 

. Grilli & Horrillo (1997h. 

Grilli et al. (|200lh. Fochesatol (|2004h. Fochesato & Diad (200611. ] 

^ochesato et all (2007hb 


Its kinematical ac curacy has been validated against the analytical solutions for infinites¬ 
imal sine waves in PhillipsI (1977). 


WSIM uses a boundary element method (BEM) to compute field variables. The 16- 
node quadrilateral elements provide global third order precision, and high-order tan¬ 
gential derivatives needed for the time discretisation are computed in a local 25-node 
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(a) (b) 


Figure 1: Evolution details of the initiation of crest roll-over at the leading crest of a 
2D C3N5 deep water wave packet undergoing weak breaking. The sequence of local crest 
surface profiles is shown in panel (b), starting after the appearance of the vertical tangent 
on the forward face. The crest tip moves from left to right. The small solid circles indicate 
the computational nodes. Details of the splined curves are given in section 


quadrilateral element curvilinear coordinate system giving fourth order precision. A fast 
multipole algorithm is used to invert the BEM problem. 

3D simulations, with x as the main direction of propagation, y the transverse horizontal 
direction and z the vertical direction, were run using mainly 16 nodes per wavelength. 
The insensitivity of the results to the resolution was established by a subset of runs using 
32 nodes per wavelength. The number of nodes in the y—direction or in the z—direction 
was adjusted to keep the boundary element aspect ratio close to 1. Figure [Tal shows a 
typical elevation profile of a gently breaking 2D C3N5 deep water wave packet just after 
exceedance of our breaking onset threshold. This case is discussed in detail in section !?^ 
and shown in figures Ffbl and iTdl Figure [lb] shows the sequence of elevation profiles of this 
overturning wave crest at incremental time steps. It highlights the ensuing development 
of the crest tip jet following the occurrence of a vertical tangent on the forward face, 
confirming the overturning capability of WSIM. 

Figure m shows a typical simulation using 16 nodes per wavelength. It illustrates the 
breaking initiation of the crest of a 3D converging deep water chirped wave packet with 
5 waves in the temporal group. Even though each cell is represented by a flat quad¬ 
rangle in the visualisations in figure [51 bi-cubic and fourth-order basis polynomials were 
used, respectively, for the physical variables and for the geometry, providing second-order 
curvature discretisation. 

Each breaking case we investigated conforms to this generic systematic progression 
(exceeding the breaking threshold, subsequent vertical tangent on forward face of crest 
and formation of initial crest tip jet). The numerics handles this smoothly, with the code 
stopping when the boundary elements at the crest tip jet become enmeshed. However, 
this occurs well beyond the exceedance of the breaking onset threshold and has no impact 
in determining this threshold. This is described in section 13.41 and in greater detail in 
appendix IA.3I 

The mixed Eulerian-Lagrangian numerical scheme and the high nonlinearity of the 
waves make the free surface mesh prone to distortion by the Lagrangian drift current. 
Extreme care was taken so that even at maximum recurrence, only a moderate Lagrangian 
drift was produced and the mesh did not deform significantly. An important benefit of 
the local Lagrangian drift is the clustering of the nodes around sharp crests. 
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Figure 2: Numerical wave tank simulation showing the initiation of breaking at the leading 
crest of a 3D converging C3N5 deep water wave packet as it propagates to the right. 
Panel (a) shows the computational domain while panel (b) shows a zoomed-in view of 
the overturning crest. The simulation used 16 nodes per wavelength, a paddle amplitude 
of Ap/Xp = 0.2067 and a linear focal distance of x/Xp = 5. The rectangular cells are 
artefacts of the visualisation. 


3.3. Scaling parameters in the numerical simulations 
The numerical simulations use the following deep water (DW) non-dimensional param¬ 
eter scalings: reference length scale is the wavelength Xp at the wavemaker (with cor¬ 
responding wavenumber fcp); reference time scale is based on the baseline frequency uip 
of the paddle. The deep water dispersion relation imposes a gravitational acceleration 
gow = 27r and a reference linear phase speed cdw = 1- Results have been produced for 
depth to wavelength ratios {d/Xp) belonging to the interval [0.2, Ij. 

3.4. Numerical convergence 

This research aims to identify a criterion that can robustly predict whether growing 2D 
or 3D carrier waves in evolving wave groups will attain their maximum steepness without 
breaking (recurrent waves) or proceed to break, with overturning crests. To achieve this 
aim, the weakest form of breaking, marginal breaking, is computed and compared with 
the corresponding maximum recurrent case. We carried out a detailed sensitivity study 
that demonstrates convergence of the NWT model even for such low-intensity breaking. 
Sensitivity tests were performed for one 2D ensemble, C3N5, using the two different 
resolutions described above (16 and 32 nodes per wavelength). We confirmed that both 
resolutions share the initial stage of breaking onset, conhrmed by the occurrence of 
a vertical forward face segment and a multiple-valued free surface. While some minor 
differences were seen between the different resolutions, the numerical convergence of 
the computed breaking parameter is discussed in detail in appendix IA.31 where it is 
established that the resolution of 16 nodes per wavelength suffices to robustly quantify 
the breaking onset threshold. 


4. Analysis of simulation data using previous breaking criteria 

The recent review paper of IPerlin et all ( 2013ll provides a detailed analysis of the 
different classes of proposed breaking criteria. Shortcomings have been identified in each 
of the criteria proposed to date. In the present study, several of these breaking criteria 
were also investigated using our simulation data and their validity evaluated. 
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Breaking 

Non-breaking 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 


Figure 3: Distribution of local crest maximum steepness Sc partitioned according to 
breaking or recurrent events, for all the C3 wave packets investigated in this study. Note 
the extensive Sc range shared by non-breaking and breaking waves. 


4.1. Geometrical criteria 

As reported in Perlin et all ( 2013h . geometric threshold criteria were not found to be 
robust (in the generic sense) in previous observational studies. Figure [3] provides an 
overview of the performance of the steepness criterion for breaking onset for the present 
data set. It shows the distribution of the maximum local crest steepness Sc = TrajXc of 
each of the recurrent or breaking wave packets. Here a is the crest amplitude and Ac is 
the horizontal exten d, measured b e tween the wave profile zero up and down crossings 
that span the crest ( Banner et all (2014)). Our results confirm previous findings that 
breaking onset cannot be discriminated by the steepness criterion since some breaking 
crest cases have a significantly lower Sc than recurrent waves. 

4.2. Kinematic criteria 

As described in Perlin et al. (1201 3 ), various k inema tic b reaking criteria have be e n pro - 
posed. We note that Stansell fc MacFarlanj ( 2002[l and Kurnia fc van GroesenI (l2014h 
mention a variety of purely kinematic criteria (ratio of fluid speed to wave speed larger 
than 0.7, 0.8 and 1) which are embraced by our findings in the present paper. However, 
there is no guidance as to the u niversality of the detection of breaking onset. 

Recently. [Banner et al. ( 2014l l investigated the influence of unsteadiness in dispersive 
wave packets, based on results from numerical simulations and complementary laboratory 
and ocean tower measurements. Their study highlighted the existence of a significant 
kinematic/geometric phenomenon attributable to the extra degrees of freedom in such 
unsteady wave packets. This results in an additional generic oscillatory crest/trough 
leaning mode characterising the carrier wave evolution. For focusing deep water wave 
packets, this manifests as a systematic crest speed slowdown of approximately 20% of 
the linear phase velocity, reconciling why their breaking crests are observed to have 
initial speeds typically 20% lower than the corresponding linear phase s peed for that 
wavelength. The source of the cr est s lowdown mech anism is investigated in [Banner et al 


{ 2 OI 4 I I , [Barthelemy et al\ ( 2015 ) and Fedeld ( 2014[ l . Most significantly, it determines the 
underlying crest motion from which the breaking onset initiates, which is a key element 
in our new breaking framework. 


4.3. Dynamical criteria 

Dynamical criteria link the physics of breaking onset to the energy fluxes associated with 
the underlying unsteady wave group structure. Conceptually, the rate of convergence 
of intra-group energy flux exceeds a local stability level at a particular crest, which 
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triggers this crest to break. The highly nonlinear, spatially non-uniform and unsteady 
nature of the flow field makes rigorous analysis difficult. One of the key results in section [T] 
describing dynamical breaking criteria was the proposition tha t breaking onset is i nitiat ed 
when the crest particle speed exceeds the linear group speed ( Tulin fc Landrini I 2000ll b 
We were able to investigate this criterion for chirped nonlinear wave packets for a range of 
depth conditions, including 3D cases. However, based on our findings in section [5] below, 
our results do not support their proposed criterion. 

We also revi sited the wave group-r elated energy flux approach invest igated in the mod¬ 


elling study of ISong fc: Banneil (|2002ll and validated observationally in [Banner fc Peirson 


( 2003 ) . This approach examined the possibility of a generic local parametric energy con- 
ve rgence rate for 2D wave groups following the wave group maximum (see section 3.3 
Perlin et al\ (12013)). F urther support was recently reported for 2D breaking onset by 


Derakhti fc KirbX i 201(t) (section [T]) . However, two factors motivated our present search 
for an alternative breaking criterion based on an energy convergence rate threshold. 

Firstly, our present study found that 2D and 3D break ing onset behaviour di d not 
closely match the 2D threshold growth rate proposed bv ISong fc BannCT ( 2002 1. For 
both 2D and 3 D cases in our study, t he computed local carrier wavenumber used in 
constructing the Song fc Banneil ( 2002ll growth rate did not increase monotonically as 
the wave steepened, contrary to the analysis of Song fc Banner ( 2002h . This departure 
was also observed experimentally in section 6.3 of AllisI ( 2013 1. As a resul t, the diagnostic 
growt h rate trajectory departed si gnificantly from the re sults reported in ISong fc Bann^ 
( 2 OO 2 II . Secondly, the approach of ISong fc Banner ( 2002 1 requires tracking, for any given 
crest, the space-time locus of its maximum elevation for at least 2 cycles prior to its 
reaching its ultimate maximum (either the recurrence maximum or breaking onset). 
Aside from its measurement complexity, this approach becomes tenuous for cases when 
more rapid approach to breaking onset occurs with fewer than 3 growth cycles. These 
factors underpinned our systematic search for a less restrictive energy flux-based breaking 
criterion. 


5. New breaking criterion based on the local energy flux velocity 

Conceptually, the onset of breaking may be regarded as the inability of the waveform 
to accommodate a local wave energy flux which exceeds that in the corresponding maxi¬ 
mum recurrent case. It is observed that breaking of the dominant waves typically occurs 
at the crest of the tallest dominant wave within a group, showing the preferred crest 
localisation of the phenomen on. This is consistent with the open ocean observations of 
Holthuiisen fc Herberd ( 1986h . Excess local wave energy flux can arise from a variety of 
sources, such as intra-wave energy exchanges, wind-wave exchange, geometrical and tem¬ 
poral 3D wave focusing, wave-current interactions, among others. For the present focus 
on unforced water wave groups, we hypothesise that the same breaking onset physics 
could apply whether the wave group is evolving in deep water or in intermediate water 
depth. The shortcomings of the various criteria described above led us to focus on the role 
of excess wave energy flux as the underpinning element of a gene r ic bre aking criterion. 
The wave energy flux vector is defined in section (2.3) of Philliosl (197^, together with 
a discussion of its local conservation equation (2.3.2). 
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5.1. Energy flux considerations in nonlinear wave groups 

The mechanical energy balance equation relates the local rate of change of the energy 
density E 

E = pg {z - zo) + ^p\\u\\^ 
to the divergence of the local energy flux F 

F = u(^{p- po) + pg{z- zo) + ip||uf 

where p is the pressure, po is the ambient pressure above the surface (taken as zero 
without loss of generality), ||u|| is the fluid speed, g is gravitational acceleration, z is the 

vertical coordinate and zq the datum. _ 

With the above dehnitions, equation (3.6.14) in IPhillio^ (1973) shows the conservation 
law for the depth-integrated energy. Based on (3.6.14), we performed an analysis of the 
crest behaviour of the depth-integrated energy density, depth-integrated energy flux and 
its gradient at maximal focusing for representative nonlinear wave packets. We investi¬ 
gated examples of chirped packets with different numbers of carrier waves. Our aim was 
to determine whether breaking onset provides a distinctive signature within the depth- 
integrated energy context. We concluded that this depth-integrated approach obscures 
this apparently highly localised crest instability. To detect the transition to breaking 
above the background wave energy, it became evident that a local energy flux analysis 
in the neighbourhood of the crest was needed. This is described in the following section. 


5.2. Breaking criterion based on the local energy flux velocity 

For the present p urposes, in an inertial frame of reference, the local energy density 
conservation law (I Phillina 1977 1. equation (2.3.2)) takes the Eulerian form: 


dE 


+ V-{u (E+p)) =0 


(5.1) 


From equation (15.Ill , in an inertial frame of reference, the energy flux velocity is seen 
to be the fluid velocity u. However, to gain a refined understanding of the energy flux 
to the tallest crest in an evolving nonlinear wave packet, the corresponding conservation 
e quation for a control volume moving with the (unsteady) crest velocity c has the form 
1( Tulin 2007 1. equation (1.2)) 


DcE 

Dt 


-I- V • ((u — c) E + up) =0 


(5.2) 


where 


DcE dE 


Dt dt 
free surface, equation (15.21) reduces to 


-I- c • VF is the rate of change following the crest. Along the unforced 


DcE 

Dt 


-f V- ((u-c) E) =0 


(5.3) 


which shows that (u — c) is the relevant flux velocity transporting energy to the growing 
crest. 

These theoretical aspects taken in context with allied observational and computational 
considerations, motivated our investigation of the behaviour of the local energy flux F in 
relation to the local energy density E in the neighbourhood of the wave crest in the fixed 
frame of reference. Here equation dsm can be used, along with the incompressibility 
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condition V • u = 0, to quantify the ratio F/E as follows: 

F/E = u{E + p) /E (5.4) 

On the free surface, F/E reduces to the surface fluid velocity u. Given the intrinsic 
relevance of the crest velocity c highlighted above, we adopt c as the natural normalising 
velocity for the flux ratio in equation ()5.4I) and introduce a breaking onset threshold 
parameter B as 


B = F/(i:;||c||) 


(5.5) 


The above analysis indicates the crest speed is a key variable in the breaking on- 
set threshold param eter B, as mentioned above in section 14.21 The recent study by 
Banner et all { 2014h and Barthelemv et al. { 2015h provides new insights into crest speed 


behaviour in the unsteady evolution of nonlinear dispersive water wave packets. Specif¬ 
ically, every wave in an unsteady dispersive wave group experiences a dynamic leaning 
cycle. Crests and troughs enter at the rear of the wave group hrst leaning forward, then 
transitioning through symmetry and subsequently leaning backward as they propagate 
towards the front of the wave group. In deep water, this leaning cycle creates a sys¬ 
tematic crest slowdown to about 80% of the linear phase velocity cq, with only a weak 
dependence on the (local) steepness of the crest or trough. Here Cg is the linear phase 
speed corresponding to the peak frequency ujq of the local frequency spectrum of the 
wave packet, assuming the linear dispersion relation. However, in shallower water depth, 
the waves are less dispersive and the crest slowdown effect reduces, with crest speeds 
approaching the phase speed. As will be seen below, using the crest speed c underpins 
our key finding of a robust breaking onset threshold. 

During the crest life cycle, the local energy flux speed becomes maximal near the crest 
point, as will be shown below. Since F and c are vectors, we have two convenient choices 
to construct their ratio: either to project along the wave propagation direction (taken 
here as the x—direction) or to use norms. This leads to the two following dimensionless 
quantities: 


Bx = F,/(Ac,) or H= ||F||/(A||c||). (5.6) 

Since it is found that there is no difference between the two ratios when the crest reaches 
its maximum (the vertical components of the flux and of the crest speed vanish), we 
only explore the validity of Bx as a breaking threshold parameter that embraces wave 
kinematics and energetics. As discussed in detail below, the distributions of E, F and 
B provide key insights into the surface and subsurface manifestation of breaking onset. 
We note that the bottom depth topography does not play an explicit role in the above 
discussion of our breaking onset threshold. 


5.3. Breaking onset threshold and its surface signature 

Here we report results investigating our breaking onset criterion in terms of a parametric 
threshold band, as described in detail in section [5Tl This band is bounded below by 
a maximal non-breaking condition, and above by a marginal breaking condition, for 
individual waves in the local carrier wave system. Fortuitously, the threshold parameter 
specified in equation (15.61) has a simple signature along the free surface. 

As explained above, the pressure p at the free surface is assumed constant and is taken 
as zero without loss of generality. The energy flux velocity then reduces to the fluid 
particle velocity and our proposed breaking criterion for Bx , based on excess energy flux 
of the marginal breaking case over the corresponding maximal recurrent case, remarkably 
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(a) Paddle amplitude Ap/Xp = 0.232 (b) Paddle amplitude Ap/\p = 0.237 

Figure 4: Breaking criterion B^ as a function of time for recurrent C3N9 2D chirped 
wave packets. Each trajectory curve shows the time evolution of the breaking parameter 
Bx following each carrier wave crest maximum during the packet evolution as it grows, 
attains its maximum steepness without breaking and then decays. 


reduces to a kinematic criterion at the free surface: 

Bx = Fx/Ecx = Ux/Cx > threshold (5.7) 

Accordingly, we performed a suite of numerical simulations to investigate the behaviour 
of Bx at the free surface, as defined in equation 15.71 As described in detail below, the 
onset of breaking was found to occur when 0.85 < Bx < 0.86 for our entire ensemble of 
numerical experiments addressing 2D and 3D, deep-water and intermediate-depth cases 
over flat bottom topography. In addition to this remarkable result, it is noteworthy that 
after factoring out the reduced crest speed in deep water (rs 0.8co), this corresponds to 
u/cq ~ 0.68, which is well below the often-quoted classical kinematic breaking criterion 
u/cq > 1 . 

In regard to 3D breaking validation, energy flux is a vector quantity and our criterion 
should be able to accommodate the additional lateral energy flux in 3D converging cases. 
While our 3D numerical simulations are limited to a single set (C3N5), it is a represen¬ 
tative case with very strong lateral convergence, with a focal distance of approximately 
five carrier wavelengths. 

While the surface-based Bx criterion provides the most convenient operational breaking 
onset threshold criterion, the associated subsurface distribution of Bx , defined in equation 
15.61 was also investigated and representative results are reported below in section 15.41 
The discussion now addresses the generic behaviour of Bx at the free surface. 

Representative behaviour of Bx is shown in figures ISl I4bl [5al and I5bl Figures l4bl and 
[5a1 show the time evolution of the breaking parameter Bx following each crest maximum 
during the group propagation. The two cases illustrate near-maximum recurrence and 
marginally breaking 2D deep-water packets. The hatched zone is the identified threshold 
level of 0.85 — 0.86 determined by our entire ensemble of simulations, above which all 
crests proceeded to break. In the near-maximum recurrence case (figure I4bp , the trajec¬ 
tories never cross the hatched zone, in contrast with figure [5al which shows the marginal 
breaking case of the same wave group class, where the trajectory of Bx clearly crosses 
the threshold. 

In these examples, the computed spline passing through all the local crest maxima 
also shows their longer-term evolution trajectory. In the breaking case, this spline crosses 
the hatched breaking threshold zone before breaking subsequently initiates, providing 
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Figure 5: Breaking parameter as a function of time for breaking C3N9 2D chirped 
wave packets. Each trajectory curve show the time evolution, up to breaking onset, of 
the breaking parameter following each carrier wave crest maximum during the packet 
evolution. The trajectory of the breaking crest clearly crosses the proposed breaking 
onset threshold [0.85,0.86]. 


advance warning of up to half a carrier wave period. This behaviour is representative of 
all investigated breaking cases. 

Our key results for the proposed breaking onset threshold are based on B^^ derived from 
the ensemble of systematic numerical simulation cases shown in table [1] For each generic 
packet type, the results are ordered according to increasing paddle amplitude, with the 
transition from maximal non-breaking to marginal breaking highlighted. The C3N5 re¬ 
sults at the top of the table confirm the insensitivity of the results to the resolution, as 
discussed in detail in section IA.3I of the appendix. 

Figure [HI the key figure in this paper, provides a comprehensive summary of the per¬ 
formance of the breaking onset threshold for the present data set. The breaking onset 
parameter B^ for every crest in each packet evolution is plotted against the correspond¬ 
ing crest steepness Sc- The figure shows recurrent cases at their maximum height and 
breaking crests at their onset. Both 2D and 3D, deep and intermediate water depth cases 
are included. The vertical hatched zone on the right represents the classical Stokes local 
steepness limit expressed in terms of Sc rather than ak. The horizontal hatched zone 
[0.85 < Bx < 0.86] is the breaking threshold determined from our ensemble of numeri¬ 
cal simulations. This figure highlights two significant findings. The major finding is the 
discovery of a clear separation between recurrent and breaking crests. For recurrence 
cases, our proposed breaking onset parameter Bx is always less than 0.85, below which 
crests were never found to break. However, Bx is always greater than 0.86 when breaking 
occurs. This identifies the breaking onset threshold zone for Bx as [0.85,0.86]. 

The other finding concerns the relationship between crest steepness at maximum wave 
height and breaking onset. Once again it is seen that crest steepness is clearly not a 
threshold variable that is able to discriminate between breaking and recurrence. This is 
evident as the local steepness for breaking crests can be lower than the local steepness 
of recurrent evolution cases for a given depth (or kd). We note that a straight line can 
be drawn between (0,0) and (•S’c^^,., 0.855). All the symbols for the individual crest are 
closely aligned at low steepness and then a departure from this trend occurs due to 
nonlinearities. The Bx parameter then grows faster than the steepness. This transition 
occurs at higher steepness levels in intermediate water depth than in deep water. The 
limiting case for shallow water conditions {kd 0) seems to converge towards the deep 
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S, 


Figure 6: Breaking parameter B^ plotted against local steepness Sc- Each point was ob¬ 
tained by tracking the maximum B^ for every crest in each packet evolution documented 
in table[TJ The horizontal hatched zone at 0.85 < B^ < 0.86 is the threshold which segre¬ 
gates breaking from non-breaking cases; the vertical hatched zone Sc > 0.72 is the deep 
water Stokes limit. Hollow symbols represent recurrent crests and solid symbols repre¬ 
sent breaking crests. The wave group families are labelled as follows:- circles: C3N5 2D 
deep water; downward-pointing triangles: C3N5 3D deep water; upward-pointing trian¬ 
gles: C3N7 2D deep water; leftward-pointing triangles: C3N9 2D deep water; rightward¬ 
pointing triangles: C3N9 2D 0.2 < d/Xp < 1; squares: C3N9 2D d/Xp = 0.2; pentagons: 
C3N9 3D 0.2 < d/Xp < 1. 


water Stokes limiting steepness. Table [T] summarises the cases processed for this study, 
with their properties. 


5.4. Subsurface energy flux considerations 

The breaking parameter B is based on the local wave energy flux ||F|| normalised by 
the local wave energy density E and the wave crest speed ||c||, defined at the maximum 
elevation of the crest. The energy density and the energy flux are both well-behaved 
field quantities and reach their maximum at the crest point maximum. Because these 
quantities include the potential energy, their value is defined relative to an arbitrary 
constant, the datum zq. The total energy E = pg{z — zo) + puf /i may vanish locally, 
depending on the choice of the datum zq. Accordingly, to avoid spurious singularities of 
B and their effect on the B distribution arising within the flow domain, the datum level 
zo needs to be chosen outside the flow domain. We chose the datum as twice the depth 
of the flow domain, and verified that the subsurface B distributions were insensitive to 
Zq lower than this. 

B exists both on the surface and in the interior of the wave domain. It can be computed 
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Name 

d/ \p 

Nodes / \p 

Ap/Ap 

Sc Max 

Breaking 

Bx 

Figure 

C3N5A0.05 

1 

16 

0.02521 

0.03 

No 

0.012 


C3N5A0.3 

1 

16 

0.1513 

0.206 

No 

0.237 


C3N5A0.514 

1 

16 

0.2592 

0.519 

No 

0.575 


C3N5A0.516 

1 

16 

0.2602 

0.527 

No 

0.606 


C3N5A0.518 

1 

16 

0.2612 

0.533 

No 

0.761 


C3N5A0.51852 

1 

16 

0.26133 

0.5335 

No 

0.831 

[IQ] 

C3N5A0.51856 

1 

16 

0.26136 

0.5338 

Yes 

0.869 

US] 

C3N5A0.519 

1 

16 

0.2617 

0.534 

Yes 

0.8717 


C3N5A0.53 

1 

16 

0.2673 

0.56 

Yes 

1.03 


C3N5A0.56 

1 

16 

0.2824 

0.563 

Yes 

1.059 


C3N5A0.508 

1 

32 

0.2562 

0.498 

No 

0.631 


C3N5A0.511 

1 

32 

0.2577 

0.509 

No 

0.841 

[IS] 

C3N5A0.514 

1 

32 

0.2592 

0.520 

Yes 

0.860 

□US] 

C3N5A0.516 

1 

32 

0.2602 

0.545 

Yes 

1.015 


C3N5A0.518 

1 

32 

0.2612 

0.500 

Yes 

1.128 


C3N5A0.519 

1 

32 

0.2617 

0.466 

Yes 

1.059 


C3N5A0.32X10 

1 

16 

0.1614 

0.420 

No 

0.524 


C3N5A0.33X10 

1 

16 

0.1664 

0.474 

No 

0.815 


C3N5A0.34X10 

1 

16 

0.1715 

0.588 

Yes 

0.888 


C3N5A0.35X10 

1 

16 

0.1765 

0.603 

Yes 

1.163 


C3N5A0.36X10 

1 

16 

0.1815 

0.514 

Yes 

0.910 

m 

C3N7A0.41 

1 

16 

0.2067 

0.344 

No 

0.429 


C3N7A0.42 

1 

16 

0.2118 

0.358 

No 

0.452 


C3N7A0.43 

1 

16 

0.2168 

0.374 

No 

0.480 


C3N7A0.44 

1 

16 

0.2219 

0.391 

No 

0.509 


C3N7A0.45 

1 

16 

0.2269 

0.406 

No 

0.547 


C3N7A0.46 

1 

16 

0.232 

0.428 

No 

0.588 


C3N7A0.47 

1 

16 

0.237 

0.448 

No 

0.654 


C3N7A0.48 

1 

16 

0.242 

0.484 

No 

0.737 


C3N7A0.49 

1 

16 

0.2471 

0.513 

Yes 

0.944 


C3N7A0.50 

1 

16 

0.2521 

0.523 

Yes 

0.964 


C3N9A0.42 

1 

16 

0.2118 

0.360 

No 

0.477 


C3N9A0.43 

1 

16 

0.2168 

0.375 

No 

0.510 


C3N9A0.44 

1 

16 

0.2219 

0.392 

No 

0.550 


C3N9A0.45 

1 

16 

0.2269 

0.413 

No 

0.601 


C3N9A0.46 

1 

16 

0.232 

0.437 

No 

0.667 

I4al 

C3N9A0.47 

1 

16 

0.237 

0.464 

No 

0.788 

[40 

C3N9A0.48 

1 

16 

0.242 

0.496 

Yes 

0.952 

[5al 

C3N9A0.49 

1 

16 

0.2471 

0.433 

Yes 

0.977 

[50 

D0.4C3N9A1.05 

0.2 

16 

0.529 

0.616 

No 

0.812 


D0.4C3N9A1.07 

0.2 

16 

0.54 

0.682 

Yes 

1.001 


D0.5C3N9A0.95 

0.25 

16 

0.48 

0.595 

No 

0.730 


D0.5C3N9A0.97 

0.25 

16 

0.489 

0.667 

Yes 

0.944 


DO.5C3N9A0.55X10 

0.25 

16 

0.277 

0.642 

Yes 

0.997 


D0.75C3N9A0.77 

0.375 

16 

0.388 

0.552 

No 

0.713 


D0.75C3N9A0.79 

0.375 

16 

0.398 

0.602 

Yes 

0.875 


D1C3N9A0.65 

0.5 

16 

0.33 

0.541 

No 

0.794 


D1C3N9A0.67 

0.5 

16 

0.338 

0.590 

Yes 

1.060 



Table 1: Maximum steepness Sc and maximum for each simulated wave group within 
the several ensembles listed. In each group name, the first two N[I] in the group name 
denotes the number of waves [I] in the temporal packet, Ap/Xp is the paddle amplitude 
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for maximum recurrent crests and for crests evolving to breaking, up to the point of 
breaking onset. Figure [7] and figure [8] highlight key properties of interest, respectively, 
of 2D and 3D waves evolving towards breaking. Each of these figures represents two 
snapshots (A and B) at different stages of development of the wave field. The first is a 
colour-coded representation showing contours of the a;—component of the breaking onset 
parameter (Ba:) on the crest surface (figures [Tal I7bl I5al and I8bl) . The other plot is a 
vertical plane slice on the symmetry axis of the crest, taken along the blue line shown in 
the surface plots (figures [7^ [Tdl |5^ andlSdl). 

These plots highlight the localisation of the B^ maximum at the crest point of the 
tallest wave, which is similar for the 2D and 3D cases. Of further interest is the spatial 
extent of the region involved in a breaking onset event. Based on the B^ parameter 
threshold, the horizontal extent of the zone where Bx values exceed 50% of the breaking 
onset threshold. 

Observational validation of our new breaking criterion is challenging, due to its spa¬ 
tial localisation within a rapidly-evolving, compact zone. Its horizontal extent where 
Bx exceeds the breaking threshold is only about 3% of the wavelength. Since the non- 
dimensional wave crest speed is close to 0.85 at breaking, a fixed probe will reside in 
this zone for only about 0.04 non-dimensional time units. Accurate determination of the 
rapidly-evolving crest speed and its associated surface fluid speed is a demanding mea¬ 
surement. Careful physical experiments are needed to valid ate these new findin gs. Such 
studies have been undertak en and the results r eported in Saket et~^ ( 2017ah for 2D 
waves in deep water and in ISaket et all (20176) for intermediate water depths. These 
studies encompassed different classes of nonlinear deep water and intermediate depth 
wave groups, for a range of group bandwidths. The influence of wind forcing was also in¬ 
vestigated for deep water cases. Their results show agreement with our proposed breaking 
onset threshold criterion to within 2.5%, which is close to the experimental error bounds. 

It is seen that the spatial regions where Bx becomes appreciable in 2D and 3D breaking 
waves have approximately the same vertical extent, but the surface distribution of Bx 
differs considerably. For 3D converging waves, the maximal values are positioned on side 
lobes off the symmetry axis. As the wave steepens, these lobes converge towards the 
symmetry axis, with Bx maximising on the symmetry axis where the wave breaks. In 
contrast, 2D breaking waves have an almost constant spanwise distribution of Bx values, 
with only minor variations found where the loci of maximum Bx values are on either side 
of the symmetry axis. 


5.5. Possible influence of surface tension on the energy fluxes 

The above investigation was carried out neglecting surface tension effects, which may be¬ 
come significant in zones of higher surface curvature. We assessed the validity of this as¬ 
sumption at the crest point of a 2D C3N5 wave transitioning through the Bx = [0.85, 0.86] 
breaking threshold with a wavelength of order 1 m. We found the radius of curvature 
/ 

i?=(I-|-(<C') J /(("siO.Ol, where z = ( specifies the free surface. The corresponding 
non-dimensional water-side pressure increment 6p = a/pR Ri 0.008 is to be compared 
in the Bernoulli equation to the non-dimensional kinetic energy KE = u^/2 m 0.3 and 
non-dimensional potential energy PE = gC ^ 0.58 for the datum taken at Zq = 0. 

Thus the maximum water-side pressure increment associated with surface tension is 
estimated to be about 2 orders of magnitude lower than the energies, and therefore makes 
negligible difference to the energy fluxes. Further, we note that crest point curvature of 
our numerical wave profiles as they transition the breaking onset threshold, which is well 
in advance of visual crest overturning, are only about three times the observed crest 
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curvatur e of O f 1 m) wavelengt h laboratory waves just before they proceed to overturn 
(e.g. see lOiao fc DuncanI (j200l[ l. Figure 5(a)). 

For water waves with even shorter wavelengths, there is the potential for our proposed 
breaking onset threshold to be modified by the effects of surface tension. This will only 
be resolved through detailed investigation. 


6. Comparison with previously proposed breaking onset criteria 

Our results do not support the criterion described bv iTulin fc Landrini I 2000[l which 
links breaking onset to the crest fluid speed exceeding the linear group speed of the wave 
packet. 

Based on our results, for deep water waves their threshold is equivalent to Bx = 0.62, 
which signals breaking onset at a considerably lower value than our Bx threshold. How¬ 
ever, in section 15.41 we highlighted the very close agreement (w ithin 2.5%) between 
our Bx threshold value and the measured threshold reported by ISaket et al. (2017a) 
and Saket et al. ( 2017&I1 1 for modulationally breaki ng deep water and intermediate wa¬ 
ter depth wave packets. Further, as reproduced in Saket et~^ ( 2017all . our very steep 
deep water recurrent waves ro utinely exceed the proposed group velocity-based thresh¬ 
old of Tulin fc Landrinil ( 2000ll without breaking. For intermediate water depth waves, 
as the linear group spee d approaches the linear phase speed, closer agreement with the 
Tulin fc Landrini (I 2 OOOII result might result for a specific water depth/wavelength ratio. 


Perlin et al\ ( 201.‘lt l cite two studies that carefully c ompare the highest water speeds 


with coincident wave speeds. Just prior to overturning, Perlin et al. ( 1996h found a ratio 
of maximum water speed to linear wave speed of 0.74 for a single ca se of deep water 
plunging breaking. For intermediate water depths, Chang fc Liul ( 1998ll found a ratio of 
m aximum water speed to linear wave speed of 0.86 prior to breaking. However, as shown 
bv [Banner et all ( 2014h . there is a systematic slowdown of the crest point as it transitions 
through its local maximum and this crest slowdown was not taken into account by these 
investigators. After the crest slowdown is taken into account, their measurements are 
consistent with our findings, as both investigations reported water speeds exceeding 0.85 
of the crest speed at breaking onset, with subsequent water speeds exceeding wave speeds 
d uring the overturnin g breaking process. 

Saket et al. ( 2017ah summarise results from a suite of other investigations that com¬ 


pared crest water speed with the cre st speed just prior to breaking onset. The most im¬ 
portant of these studies is t he work of Stansell fc MacFarlana ( 2002 1. whose findings were 
found to be consistent with lSaket et aZ.I ~ 2017all . once plausible and appropriate correc¬ 
tions w ere made to the near-surface velocity structure reported by Stansell fc MacFarlanel 

(I 2 OO 2 II . 


7. Discussion and conclusions 

A new breaking onset threshold criterion based on energy flux considerations has been 
developed for water waves and its applicability investigated for 2D and 3D chirped fo¬ 
cusing wave groups in deep and intermediate water depths. While the energy flux that 
initiates breaking arises from energy focusing within the whole wave group, we found 
that the initial instability occurs within a very compact region, about 3 percent of the 
wavelength in horizontal extent, centred on the maximum wave crest. 

We found that depth-integrated quantities were not able to detect the signature of 
breaking generically in the phase-resolved wave motion. A more localised examination of 
the flow near the free surface was necessary to detect the initiation of breaking and to find 
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a consistent deterministic indicator. Our new breaking criterion is based on the strength 
of the local energy flux relative to the local energy density, normalised by the local crest 
speed. This non-dimensional parameter B, which reduces to B^ at the crest point (here 
X is the direction of propagation), simplifies fortuitously for zero surface pressure. On 
the free surface, it reduces to a kinematic condition for the ratio of crest fluid speed to 
crest point speed. 

For the condition of zero surface pressure forcing, our suite of numerical experiments 
using a significant range of chirped wave packets has shown that our new dynamically- 
based breaking onset parameter B^ has a generic threshold value computed to be in the 
range [0.85, 0.86], above which any local crest will undergo breaking onset. This threshold 
band is applicable everywhere in the fluid domain and has important new consequences. 

Surface projection of our criterion is needed to make comparisons with previously 
proposed kinematic criteria. When projected along the surface, our criterion predicts 
breaking onset considerably in advance of the kinematic breaking criterion, where break¬ 
ing is associated with crest fluid speed exceeding the crest point speed. For maximally 
steep non-breaking waves, the crest fluid speed cannot exceed 0.85 of the crest speed. 
When referenced to the linear wave speed Cq, after factoring in the generic crest slow¬ 
down for focusing dispersive deep water wave packets for which c sa 0.8co (see sections 
14.21 and 15.2|) . our B^ threshold predicts breaking onset for deep water waves when the 
crest fluid speed to linear phase speed ratio u/cq exceeds 0.68. In progressively shallower 
water depth conditions, the generic crest slowdown reduces (c —^ cq) and u/cq increases 
towards [0.85, 0.86]. In any event, it is evident that our breaking onset threshold band, 
which follows from energy flux considerations, occurs well before the water speed outruns 
the crest speed. 

Overall, remarkable robustness of this new breaking onset threshold was found for the 
diverse range of 2D and 3D chirped focusing wave packets on deep and intermediate depth 
flat bottom topography investigated in this study. The threshold is also con sistent with 


the find ings of available compar able experiments for wavelengths of O (Im) (|Saket et al 


( 2m7r7h and lSaket et all (j2ni7lJ ll. Further investigation is needed for shorter waves when 
surface tension effects may become important and also for other commonly encountered 
breaking wave scenarios, such as waves shoaling on inclined bottom topography, coalesc¬ 
ing wave groups and waves on uniform and shear currents. For these scenarios, since there 
is no explicit dependence on water depth, mean current or energy convergence rate in our 
energy flux considerations (section [Ol , we anticipate that our breaking onset threshold 
will be applicable to such cases. 
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(A) 



(a) Colour-coded zones of the breaking parameter on the free surface at t = 13.85. The 
blue line is the location of the slice seen in figure 


(B) 



7.5 


(b) Same as (a) at t = 14.41. The blue line is the location of the slice seen in figure Fzdl 




(c) Colour-coded cross-section within the slice (d) Colour-coded cross-section within the slice 


y = 0 shown in figure [7al 


y = 0 shown in figure [7b] 


Figure 7: Contours of the breaking criterion for the 2D C3N5A0.514 breaking wave 
case. The regions are colour-coded as follows: red: > 0.85; green: 0.7 < B^ < 0.85; 

blue: 0.6 < B^ < 0.7; yellow: 0.5 < B^ < 0.6; magenta: 0.4 < B^ < 0.5; white: B^ < 0.4. 





20 


X. Barthelemy, M.L. Banner, W.L. Peirson, F. Fedele, M. Allis, F. Dias 


(A) 



(a) Colour-coded zones of the breaking parameter on the free surface at t = 12.45. The 
blue line is the location of the slice seen in figure 15^ 


(B) 



(b) Same as (a) at t = 12.84. The blue line is the location of the slice seen in figure I8d1 




(c) Colour-coded cross-section within the slice (d) Colour-coded cross-section within the slice 


y = 0 shown in figure [Sal 


y = 0 shown in figure [8b] 


Figure 8: Same as figure [71 but for the 3D C3N5A0.36X10 breaking wave case. 
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Appendix A. Details of the numerical wave tank 

A.l. The numerical wave tank 

WSIM uses a mixed Eulerian-Lagrangian time-updating scheme for irrotational motion 
described by the velocity potential (x, t), in a Cartesian coordinate system x = (cc, y, z) 
with constant pressure at the open water surface. 2 is the vertical upward direction and 
z = 0 the still water surface (figure [9]). 

Fluid velocity is defined as u = = (u, u, ir). The continuity equation leads to the 

Laplace equation for the potential within the fluid domain D (t) lAl]). The symbols F 
and F ps are used to denote the entire domain boundary and the free surface respectively. 
The equations read: 


V'^cj) = 0, on D 

(Al) 

Dr 

— = u = V0, on Ffs 

(A 2) 

^ = -gz + - —, on Tps 

Ut 1 p 

(A3) 

dnfi = 0, on r\rFS 

(A 4) 


where r is the position vector of a fluid particle on the free surface, g the acceleration 
due to gravity, po the atmospheric pressure, p the fluid density and (= ^ + V(/>-V) 
the Lagrangian (or material) time derivative. 

The free surface (domain Fps) is described by fully-nonlinear kinematic (jA 2ll and 
dynamic (jA 3ll equations. Recent developments have implemented a 3D snake wave paddle 
at one vertical face of the domain which is described subsequently. The far face of the 
numerical wave tank from t he paddle has an absorb ing beach which damps any incident 
wave energy as described in Grilli fc Horrillol 1 1997t l. All remaining faces of the domain 
have a zero-flux boundary condition (r\ri? 5 ) (jA 4^ . The depth of the NWT domain is 
equal to the wavelength (Ap) and its size (length, width) is (12.5Ap, l/4Ap) for 2D and 
(lOAp, 7Ap) for 3D computations. 


A.2. Post-processing 

A significant challenge within this investigation was determining quantities immediately 
below the highly-curved free surface. While the surface velocity distribution is available 
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Figure 9: Schematic (not to scale) of the WSIM simulation domain and nomenclature. 


directly from the model output, reconstructing the subsurface velocity field values re¬ 
lies on the ability to correctly estimate Green’s integral over the entire domain. Three 
complementary methods are used depending on the distance between the inner point 
and the boundaries. When the distance is large enough (relative to the element size), 
the integration on the individual element is carried out by a Riemann quadrature on 
a classical Gauss-Lobatto point distribution. When the distance becom es sma l l, this 
classical quadrature method does not retain enough precision, and the Telled ( 19871) 
method is used. The Telles method consists of binary subdivisions of the integration 
space to maximize precision where the singularity begins to manifest itself. The preci- 
sio n of this method is acceptab le at moderate distance from the boundary, as described 
in Grilli fc Siibramanval ( 1994l l. However, the Telles method becomes inefficient when 


the near-boundary singularities are too strong, specifically adjacent to the highly curved 
surface of steep waves. Consequently, a third method was developed a nd implemente d 
called Projection and Angula r and Radial Tr ansformation (PART) (see Havamil (1991), 
Havami fc Matsumotol ( 1994h . |h avam ] (l200,'iM . 


The crest and trough locations are tracked for each carrier wave in the packet using 
a two-stage detection algorithm. First, extrema and the zero-crossing positions are com¬ 
puted semi-analytically from the same order spline polynomials used in the BEM 
code. These positions are then linked together between time steps to form the space-time 
characteristics of the motion, where crests, troughs and zero-crossings were detected and 
followed in the inertial frame of reference of the wavetank. The speed of each crest was 
then computed as the first derivative of these trajectories in space and time. 


A.3. Details on the convergence of the NWT for marginally breaking waves 
In this section, a more detailed analysis of maximal recurrent and marginal breaking cases 
is made to investigate possible sensitivity of the breaking onset parameter computed using 
two different resolutions: one using an average of 16 nodes per wavelength (R16) and a 
higher resolution using 32 nodes per wavelength (R32). Small differences between the 
runs are reported and discussed below, but are of minor consequence and the NWT is 
found to be sufficiently accurate at R16 resolution to objectively resolve the differences 
between maximum recurrent and marginal breaking cases. For this purpose we adopted 
the occurrence of a vertical segment on the forward face as a consistent post-breaking 
onset reference for comparing sensitivity of the paddle amplitudes, space-time locations. 
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Nodes / Xp 

Figure 10: Breaking index Bx versus average number of nodes per wavelength Xp for the 
same C3N5 case. Filled circles indicate breaking crests, while open circles are recurrent 
crests. Each point is derived from the maximum of each individually-tracked wave as in 
figures 01 [S] and IHl The hatched breaking threshold zone based on our entire ensemble of 
derived Bx values is seen to be relatively insensitive to the resolution. 

wave steepness and fluid velocities for the 2D C3N5 marginal breaking case computed 
at different resolutions. Results for our proposed breaking parameter Bx are shown in 
figure [TOl 

The numerical experiments we ran for 2D C3N5 wave packets for different resolutions 
showed some minor differences. In this instance, the paddle amplitude necessary to reach 
maximum recurrence or marginal breaking reduces slightly between R16 and R32 cases. 
The marginal breaking results for the slightly steeper R32 run produced slightly larger 
breaker amplitude, steepness and energy flux than the R16 case, consistent with the 
observed slightly shorter breaking onset epoch and corresponding fetch. However, for our 
primary goal of computing the derived breaking parameter Bx, different resolutions for 
the same C3N5 wave group have only a minor influence on the determination of the 
threshold. Figure [TO] confirms the convergence of the NWT simulations in determining 
the threshold zone for Bx- Our proposed nominal breaking threshold value Bx = 0.85 
is well-bracketed for each resolution: Bx is insensitive to the resolution, and R16 is seen 
to be sufficient for our purpose of establishing a discriminating breaking criterion for 2D 
and 3D deep and intermediate depth water focussing wave packets. 
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